/*
* Copyright 1993-2010 NVIDIA Corporation.  All rights reserved.
*
* Please refer to the NVIDIA end user license agreement (EULA) associated
* with this source code for terms and conditions that govern your use of
* this software. Any use, reproduction, disclosure, or distribution of
* this software and related documentation outside the terms of the EULA
* is strictly prohibited.
*
*/

/*  This file includes modifications to the SDK 
 *  Several volume rendering techniques were added 
 *
 *  Authors: Jim Pelton and Vanessa Gertman
 *  CS464 Fall 2011 Final Project
 *  Prof. Alark Joshi
 *  Boise State University
 */

#ifndef _VOLUMERENDER_KERNEL_CU_
#define _VOLUMERENDER_KERNEL_CU_

#include "volumeRender_param.h"

#include <cutil_inline.h>   
#include <cutil_math.h>

typedef unsigned int  uint;
typedef unsigned char uchar;
typedef unsigned char VolumeType;

typedef struct {
    float4 m[3];
} float3x4;


struct Ray {
    float3 o;	// origin
    float3 d;	// direction
};

__constant__ float3x4 c_invViewMatrix;  // inverse view matrix
__constant__  vrParams c_kparams;

texture<VolumeType, 3, cudaReadModeNormalizedFloat> tex;         
texture<float4, 3, cudaReadModeElementType> gradientsTex;
texture<float4, 1, cudaReadModeElementType> transferTex; 

cudaArray *d_volumeArray = 0;
cudaArray *d_transferFuncArray;
cudaArray *d_gradientsArray = 0;

extern "C"
void setTextureFilterMode(bool bLinearFilter)
{
    tex.filterMode = bLinearFilter ? cudaFilterModeLinear : cudaFilterModePoint;
}

 __device__ 
int intersectBox(Ray r, float3 boxmin, float3 boxmax, float *tnear, float *tfar)
{
    // compute intersection of ray with all six bbox planes
    float3 invR = make_float3(1.0f) / r.d;
    float3 tbot = invR * (boxmin - r.o);
    float3 ttop = invR * (boxmax - r.o);

    // re-order intersections to find smallest and largest on each axis
    float3 tmin = fminf(ttop, tbot);
    float3 tmax = fmaxf(ttop, tbot);

    // find the largest tmin and the smallest tmax
    float largest_tmin = fmaxf(fmaxf(tmin.x, tmin.y), fmaxf(tmin.x, tmin.z));
    float smallest_tmax = fminf(fminf(tmax.x, tmax.y), fminf(tmax.x, tmax.z));

    *tnear = largest_tmin;
    *tfar = smallest_tmax;

    return smallest_tmax > largest_tmin;
}

// transform vector by matrix (no translation)
 __device__ 
float4 mul(const float3x4 &M, const float4 &v)
{
    float4 r;
    r.x = dot(v, M.m[0]);
    r.y = dot(v, M.m[1]);
    r.z = dot(v, M.m[2]);
    r.w = 1.0f;
    return r;
}

__device__ 
float3 mul(const float3x4 &M, const float3 &v)
 {
     float3 r;
     r.x = dot(v, make_float3(M.m[0]));
     r.y = dot(v, make_float3(M.m[1]));
     r.z = dot(v, make_float3(M.m[2]));
     return r;
 }

 __device__ 
uint rgbaFloatToInt(float4 rgba)
{
    rgba.x = __saturatef(rgba.x);   // clamp to [0.0, 1.0]
    rgba.y = __saturatef(rgba.y);
    rgba.z = __saturatef(rgba.z);
    rgba.w = __saturatef(rgba.w);
    return (uint(rgba.w*255)<<24) | (uint(rgba.z*255)<<16) | (uint(rgba.y*255)<<8) | uint(rgba.x*255);
}

 /*
  *	d_depth_buffer
  * create depth information 
  */
  __device__ 
  void d_depth_buffer(float *d_dbuf, int x, int y, float sample, int maxSteps, int i)
  {
     // float depth = length(pos - eo);
      float depth = (i/maxSteps)*255.f;
      d_dbuf[y*c_kparams.imageW + x] = depth;
  }

/* Classification Computation, used in d_render_classify
 * Technique Adpated from: Display of Surfaces from Volume Data
 *                         by Marc Levoy
 * classification via gradient magnitude and scalar value of data.
 * returns alpha value for given sample.
 */
 __device__
float classify(float sample, float gradMag)
 {

     float alpha = 0.f;
     float total = 1.f;
     float *alpha_n = c_kparams.isoAlphas;
     float *iso_n = c_kparams.isoValues;
     int numvals = c_kparams.numIsoValues;
      for (int i = 0; i < numvals; i++)
      {
         if (gradMag==0 && sample==iso_n[i])
         {
             alpha = 1;
         }
         else if (    gradMag > 0
             &&  sample - c_kparams.fallOff*gradMag <= iso_n[i]
             &&  iso_n[i] <= sample + c_kparams.fallOff*gradMag)
         {
             alpha  = alpha_n[i] * (1 - (1/c_kparams.fallOff)*( fabsf( (alpha_n[i] - sample)/(gradMag) )));
         }
         total = (total * (1-alpha));
     }
     return 1-total;
 }
/*
*	initCuda(void *h_volume, float *h_gradmag, float *h_gradients, 
*            cudaExtent volumeSize)
*
*  void  *h_volume:       scalar data from the host
*  float *h_gradients:    gradients from the host
*  cudaExtent volumeSize: size of x,y,z dims of h_volume.
*/
extern "C"
void initCuda(void   *h_volume, 
              float  *h_gradients, 
              float4 *h_transFunc,
              float  *d_dbuf, 
              size_t tfuncSize, 
              cudaExtent volumeSize)
{

    /************************************************************************
    * Scalar data texture                                                                  
    ************************************************************************/
    cudaChannelFormatDesc channelDesc = cudaCreateChannelDesc<VolumeType>();
    cutilSafeCall( cudaMalloc3DArray(&d_volumeArray, &channelDesc, volumeSize) );

    cudaMemcpy3DParms copyParams = {0};
    copyParams.srcPtr   = 
        make_cudaPitchedPtr(h_volume, volumeSize.width*sizeof(VolumeType), 
        volumeSize.width, volumeSize.height);
    copyParams.dstArray = d_volumeArray;
    copyParams.extent   = volumeSize;
    copyParams.kind     = cudaMemcpyHostToDevice;
    cutilSafeCall( cudaMemcpy3D(&copyParams) );  

    tex.normalized = true;                     
    tex.filterMode = cudaFilterModeLinear;     
    tex.addressMode[0] = cudaAddressModeClamp; 
    tex.addressMode[1] = cudaAddressModeClamp;

    cutilSafeCall(cudaBindTextureToArray(tex, d_volumeArray, channelDesc));



    /************************************************************************
    * Gradient Vector texture                                                                 
    ************************************************************************/
    cudaChannelFormatDesc gradientsArray_channelDesc = cudaCreateChannelDesc<float4>();
    cutilSafeCall( cudaMalloc3DArray(&d_gradientsArray, 
        &gradientsArray_channelDesc, 
        volumeSize) );

    cudaMemcpy3DParms gradientsArray_copyParams = {0};
    gradientsArray_copyParams.srcPtr = 
        make_cudaPitchedPtr( h_gradients, volumeSize.width*sizeof(float4), 
        volumeSize.width, volumeSize.height);
    gradientsArray_copyParams.dstArray = d_gradientsArray;
    gradientsArray_copyParams.extent   = volumeSize;
    gradientsArray_copyParams.kind     = cudaMemcpyHostToDevice;
    cutilSafeCall( cudaMemcpy3D(&gradientsArray_copyParams) );  

    gradientsTex.normalized = true;
    gradientsTex.filterMode = cudaFilterModeLinear;
    gradientsTex.addressMode[0] = cudaAddressModeClamp;  
    gradientsTex.addressMode[1] = cudaAddressModeClamp;

    // bind array to 3D texture
    cutilSafeCall(cudaBindTextureToArray(gradientsTex, d_gradientsArray, gradientsArray_channelDesc));
     cudaChannelFormatDesc transfunc_channelDesc = cudaCreateChannelDesc<float4>();
     cudaArray *d_transferFuncArray;
     cutilSafeCall(cudaMallocArray( &d_transferFuncArray, &transfunc_channelDesc,tfuncSize,1));
     cutilSafeCall(cudaMemcpyToArray( d_transferFuncArray, 0,0, h_transFunc, tfuncSize*sizeof(float4), cudaMemcpyHostToDevice));

    transferTex.filterMode = cudaFilterModeLinear;
    transferTex.normalized = true;
    transferTex.addressMode[0] = cudaAddressModeClamp;

    cutilSafeCall( cudaBindTextureToArray( transferTex, d_transferFuncArray, transfunc_channelDesc));
}


extern "C" 
    void freeCudaBuffers()
{
    cutilSafeCall(cudaFreeArray(d_volumeArray));
    cutilSafeCall(cudaFreeArray(d_transferFuncArray));
    cutilSafeCall(cudaFreeArray(d_gradientsArray));
}

extern "C"
    void copyInvViewMatrix(float *invViewMatrix, size_t sizeofMatrix)
{
    cutilSafeCall( cudaMemcpyToSymbol(c_invViewMatrix, invViewMatrix, sizeofMatrix) );
}

extern "C"
void copyKParamsToDevice(struct vrParams *par, size_t sizeofParams)
{
	cutilSafeCall( cudaMemcpyToSymbol(c_kparams, par, sizeofParams));
}

extern "C"
void copyDepthBufferToHost(float *h_dest, float *d_src, size_t sz) 
{
    //fprintf(stdout, "copyin' it\n");
    cutilSafeCall( cudaMemcpy(h_dest, d_src, sz, cudaMemcpyDeviceToHost));
}

/* Distance Color Blending Computation, will be called in d_render_dcb
 * Technique Adpated from: Volume Illustration: Non-Photorealistic Rendering
 *                         of Volume Models by Ebert and Rheingans
 * Makes colors which are farther from the viewer darker, they recede in the distance
 * @param col - the sample color
 * @param i - the step through the ray
 */
__device__ float4
compute_dcb(float4 col, int i)
{
	float4 cb = make_float4(0.f, 0.f, 0.15f, 1.f);
	//compute the fractional distance through the volume and raise it to kde
	//to get fractional dist we divide i by the maxsteps, in this case we multiply i by .002 (ie i/500)	
	float dv = pow(i*0.002f, c_kparams.kde);
	col.x=(1-c_kparams.kds*dv)*col.x+(c_kparams.kds*dv*cb.x);
	col.y=(1-c_kparams.kds*dv)*col.y+(c_kparams.kds*dv*cb.y);
	col.z=(1-c_kparams.kds*dv)*col.z+(c_kparams.kds*dv*cb.z);

	return col;
}

/* Tone Shading Computation, will be called in d_render_tone
 * Technique Adpated from: Volume Illustration: Non-Photorealistic Rendering
 *                         of Volume Models by Ebert and Rheingans
 * Volume appears to be lit from a warm light, objects facing the light are yellow
 * 	  objects away from the light are blue
 * @param col - the sample color
 * @param grad - the gradient for the sample
 * @param pos - the position of the sample
 */
__device__ float4
compute_tone(float4 col, float4 grad, float3 pos)
{
	grad.w = 0;
	float4 cw = make_float4(c_kparams.kty, c_kparams.kty, 0.f, 0.f);
	float4 cc = make_float4(0.f, 0.f, c_kparams.ktb, 0.f);
   	//the arbitrary light position
	const float4 light_pos = make_float4(0.8f,-0.5f,1.f,0.f);   //TODO: make light_pos non const.

	 //calclate the light vector
	 float4 light = make_float4(
			 pos.x - light_pos.x,
			 pos.y - light_pos.y,
			 pos.z - light_pos.z, 0.f);

	 light = normalize(light);
	 float lgrad = 1.0f + dot(grad,light);
	 //linearly interpolate between the blue and yellow color, use the lgrad as the weight
	 float4 It = (lgrad/2.0f)*cw + (1.0f-(lgrad/2.0f))*cc;
	 float4 Io;
	//use a white light to calculate the diffuse component
	 if (lgrad > 0)
	 {
		 Io = make_float4(make_float3(c_kparams.ktd * lgrad), 0.f);
	 } 
	 else
	 {
		 Io = make_float4(0.f);
	 }
	//set the color with the lighting values
	 col.x = It.x+c_kparams.ktd*Io.x;
	 col.y = It.y+c_kparams.ktd*Io.y;
	 col.z = It.z+c_kparams.ktd*Io.z;

	 return col;
}
/* Lighting computation
 * Technique Adpated from: Volume Rendering by Lichtenbelt et al
 * Computes the phong lighting for the sample with a single white light source
 * @param eyedir - the direction of the eye
 * @param L - the light vector
 * @param grad3 - the gradient for the sample, without the magnitude component
 */
__device__ float3
compute_light(float3 grad3, float3 L, float3 eyedir)
{
    const float3 lightCol = make_float3(1.f, 1.f, 1.f);
    const float3 col3f_a = make_float3(1,1,1);
    const float3 col3f_d = make_float3(1,1,1);
    //const float3 col3f_s = make_float3(1,1,1);

    float3 R = (2 * dot(grad3, L))*grad3 - L;
    float RdotV = dot(R, -1 * eyedir);

    float NdotL = dot(grad3, L);
    return (
            col3f_a * c_kparams.mat_a * col3f_d +
                lightCol*(c_kparams.mat_d*NdotL) + pow(c_kparams.mat_s*RdotV, 2.f)
            );
}
/* Silhouettes Computation, will be called in d_render_sil
 * Technique Adpated from: Volume Illustration: Non-Photorealistic Rendering
 *                         of Volume Models by Ebert and Rheingans
 * Computes the silhouettes which help viewer percieve shape
 * @param eyedir - the direction of the eye
 * @param grad3 - the gradient for the sample, without the magnitude component
 */
__device__ float
compute_sil(float3 grad3, float3 eyedir)
{
    float dotprod = dot(grad3, eyedir);  
    return ( c_kparams.ksc + c_kparams.kss * pow( 1-fabs(dotprod), c_kparams.kse)) * 0.01f;
}

__device__ float
compute_boundary(float gradmag)
{
    //decrease the gradient magnitude some, to avoid
    //running out of alpha too early.
    gradmag *= 0.01;
    /*Boundry Enhancement equation*/
    return (c_kparams.kgc+c_kparams.kgs*pow(gradmag, c_kparams.kge));
}

/* Silhouettes 
 * Technique Adpated from: Volume Illustration: Non-Photorealistic Rendering
 *                         of Volume Models by Ebert and Rheingans
 * Computes the silhouettes which help viewer perceive shape
 */
__global__ void
d_render_sil(uint *d_output)
{
	const int maxSteps = 500;
	const float tstep = 0.01f;
	const float opacityThreshold = 0.95f;
	const float3 boxMin = make_float3(-1.0f, -1.0f, -1.0f);
	const float3 boxMax = make_float3(1.0f, 1.0f, 1.0f);

	uint x = blockIdx.x*blockDim.x + threadIdx.x;
	uint y = blockIdx.y*blockDim.y + threadIdx.y;
	if ((x >= c_kparams.imageW) || (y >= c_kparams.imageH)) return;

	float u = (x / (float) c_kparams.imageW)*2.0f-1.0f;
	float v = (y / (float) c_kparams.imageH)*2.0f-1.0f;

	// calculate eye ray in world space
	Ray eyeRay;
	eyeRay.o = make_float3(mul(c_invViewMatrix, make_float4(0.0f, 0.0f, 0.0f, 1.0f)));
	eyeRay.d = normalize(make_float3(u, v, -2.0f));
	eyeRay.d = mul(c_invViewMatrix, eyeRay.d);

	// find intersection with box
	float tnear, tfar;
	int hit = intersectBox(eyeRay, boxMin, boxMax, &tnear, &tfar);
	if (!hit) return;
	if (tnear < 0.0f) tnear = 0.0f;     // clamp to near plane

	// march along ray from front to back, accumulating color
	float4 sum = make_float4(0.0f);
	float t = tnear;
	float3 pos = eyeRay.o + eyeRay.d*tnear;
	float3 step = eyeRay.d*tstep;

	for(int i=0; i<maxSteps; i++) {
		float posx=  pos.x*0.5f+0.5f;
		float posy=  pos.y*0.5f+0.5f;
		float posz=  pos.z*0.5f+0.5f;

		float sample = tex3D(tex, posx, posy, posz);
		float4 gradvec = tex3D(gradientsTex, posx,posy,posz);
		float4 col = tex1D(transferTex, (sample-c_kparams.transferOffset)*c_kparams.transferScale);
		col.w *= c_kparams.density;
		
        col.w *= compute_sil(make_float3(gradvec.x,gradvec.y,gradvec.z), eyeRay.d);

		// pre-multiply alpha
		col.x *= col.w;
		col.y *= col.w;
		col.z *= col.w;

		sum = sum + col*(1.0f - sum.w); // "over" operator for front-to-back blending

		// exit early if opaque
		if (sum.w > opacityThreshold)
			break;

		t += tstep;
		if (t > tfar) break;

		pos += step;
	}
	sum *= c_kparams.brightness;

	// write output color
	d_output[y*c_kparams.imageW + x] = rgbaFloatToInt(sum);
}

/*Original volume render method from the SDK
 */
__global__ void
d_render_orig(uint *d_output)
{
	const int maxSteps = 500;
	const float tstep = 0.01f;
	const float opacityThreshold = 0.95f;
	const float3 boxMin = make_float3(-1.0f, -1.0f, -1.0f);
	const float3 boxMax = make_float3(1.0f, 1.0f, 1.0f);

	uint x = blockIdx.x*blockDim.x + threadIdx.x;
	uint y = blockIdx.y*blockDim.y + threadIdx.y;
	if ((x >= c_kparams.imageW) || (y >= c_kparams.imageH)) return;

	float u = (x / (float) c_kparams.imageW)*2.0f-1.0f;
	float v = (y / (float) c_kparams.imageH)*2.0f-1.0f;

	// calculate eye ray in world space
	Ray eyeRay;
	eyeRay.o = make_float3(mul(c_invViewMatrix, make_float4(0.0f, 0.0f, 0.0f, 1.0f)));
	eyeRay.d = normalize(make_float3(u, v, -2.0f));
	eyeRay.d = mul(c_invViewMatrix, eyeRay.d);

	// find intersection with box
	float tnear, tfar;
	int hit = intersectBox(eyeRay, boxMin, boxMax, &tnear, &tfar);
	if (!hit) return;
	if (tnear < 0.0f) tnear = 0.0f;     // clamp to near plane

	// march along ray from front to back, accumulating color
	float4 sum = make_float4(0.0f);
	float t = tnear;
	float3 pos = eyeRay.o + eyeRay.d*tnear;
	float3 step = eyeRay.d*tstep;

	for(int i=0; i<maxSteps; i++) {
		float posx=  pos.x*0.5f+0.5f;
		float posy=  pos.y*0.5f+0.5f;
		float posz=  pos.z*0.5f+0.5f;

		float sample = tex3D(tex, posx, posy, posz);
		float4 col = tex1D(transferTex, (sample-c_kparams.transferOffset)*c_kparams.transferScale);
		col.w *= c_kparams.density;

		// pre-multiply alpha
		col.x *= col.w;
		col.y *= col.w;
		col.z *= col.w;

		sum = sum + col*(1.0f - sum.w); // "over" operator for front-to-back blending

		// exit early if opaque
		if (sum.w > opacityThreshold)
			break;

		t += tstep;
		if (t > tfar) break;

		pos += step;
	}
	sum *= c_kparams.brightness;

	// write output color
	d_output[y*c_kparams.imageW + x] = rgbaFloatToInt(sum);
}

/* Classification
 * Technique Adpated from: Display of Surfaces from Volume Data
 *                         by Marc Levoy
 * classification via gradient magnitude and scalar value of data.
 */
__global__ void
d_render_classification(uint *d_output)
{
    const int maxSteps = 500;
    const float tstep = 0.01f;
    const float opacityThreshold = 0.95f;
    const float3 boxMin = make_float3(-1.0f, -1.0f, -1.0f);
    const float3 boxMax = make_float3(1.0f, 1.0f, 1.0f);

    uint x = blockIdx.x*blockDim.x + threadIdx.x;
    uint y = blockIdx.y*blockDim.y + threadIdx.y;
    
    if ((x >= c_kparams.imageW) || (y >= c_kparams.imageH)) return;

    float u = (x / (float) c_kparams.imageW)*2.0f-1.0f;
    float v = (y / (float) c_kparams.imageH)*2.0f-1.0f;

    // calculate eye ray in world space
    Ray eyeRay;
    eyeRay.o = make_float3(mul(c_invViewMatrix, make_float4(0.0f, 0.0f, 0.0f, 1.0f)));
    eyeRay.d = normalize(make_float3(u, v, -2.0f));
    eyeRay.d = mul(c_invViewMatrix, eyeRay.d);

    // find intersection with box
    float tnear, tfar;
    int hit = intersectBox(eyeRay, boxMin, boxMax, &tnear, &tfar);
    if (!hit) return;
    if (tnear < 0.0f) tnear = 0.0f;     // clamp to near plane

    // march along ray from front to back, accumulating color
    float4 sum = make_float4(0.0f);
    float t = tnear;
    float3 pos = eyeRay.o + eyeRay.d*tnear;
    float3 step = eyeRay.d*tstep;

    float3 lightPos = make_float3(c_kparams.lightx, c_kparams.lighty, c_kparams.lightz);
   

    for(int i=0; i<maxSteps; i++) {

        float posx=  pos.x*0.5f+0.5f;
        float posy=  pos.y*0.5f+0.5f;
        float posz=  pos.z*0.5f+0.5f;
        
        float sample = tex3D(tex, posx, posy, posz);
        float4 gradvec = tex3D(gradientsTex, posx, posy, posz);
        float3 grad3 = make_float3(gradvec.x, gradvec.y, gradvec.z);
        float4 col = tex1D(transferTex, (sample-c_kparams.transferOffset)*c_kparams.transferScale);

        
        float alpha = classify(sample, gradvec.w);
        if (alpha==0){
            t += tstep;
            if (t > tfar) break;

            pos += step;
            continue;
        }
        
        //col.w=alpha;
        col.w *= c_kparams.density * alpha;
        float3 col_l = compute_light(grad3, normalize(lightPos-pos), eyeRay.d);

        col.x *= col_l.x;
        col.y *= col_l.y;
        col.z *= col_l.z;

        // pre-multiply alpha
        col.x *= col.w;
        col.y *= col.w;
        col.z *= col.w;

        sum = sum + col*(1.0f - sum.w); // "over" operator for front-to-back blending

        // exit early if opaque
        if (sum.w > opacityThreshold)
            break;

        t += tstep;
        if (t > tfar) break;

        pos += step;
    }
    sum *= c_kparams.brightness;

    // write output color
    d_output[y*c_kparams.imageW + x] = rgbaFloatToInt(sum);
}

/* Boundary Enhancement 
 * Technique adapted from: Volume Illustration: Non-Photorealistic Rendering
 *                         of Volume Models by Ebert and Rheingans
 * Enhances the boundaries by increasing the opacity in areas of high gradient magnitude
 */
__global__ void
d_render_bound(uint *d_output)
{
    const int maxSteps = 500;
    const float tstep = 0.01f;
    const float opacityThreshold = 0.95f;
    const float3 boxMin = make_float3(-1.0f, -1.0f, -1.0f);
    const float3 boxMax = make_float3(1.0f, 1.0f, 1.0f);

    uint x = blockIdx.x*blockDim.x + threadIdx.x;
    uint y = blockIdx.y*blockDim.y + threadIdx.y;
    if ((x >= c_kparams.imageW) || (y >= c_kparams.imageH)) return;

    float u = (x / (float) c_kparams.imageW)*2.0f-1.0f;
    float v = (y / (float) c_kparams.imageH)*2.0f-1.0f;

    // calculate eye ray in world space
    Ray eyeRay;
    eyeRay.o = make_float3(mul(c_invViewMatrix, make_float4(0.0f, 0.0f, 0.0f, 1.0f)));
    eyeRay.d = normalize(make_float3(u, v, -2.0f));
    eyeRay.d = mul(c_invViewMatrix, eyeRay.d);

    // find intersection with box
    float tnear, tfar;
    int hit = intersectBox(eyeRay, boxMin, boxMax, &tnear, &tfar);
    if (!hit) return;
    if (tnear < 0.0f) tnear = 0.0f;     // clamp to near plane

    // march along ray from front to back, accumulating color
    float4 sum = make_float4(0.0f);
    float t = tnear;
    float3 pos = eyeRay.o + eyeRay.d*tnear;
    float3 step = eyeRay.d*tstep;
     

    for(int i=0; i<maxSteps; i++) {

        float posx=  pos.x*0.5f+0.5f;
        float posy=  pos.y*0.5f+0.5f;
        float posz=  pos.z*0.5f+0.5f;
        float sample = tex3D(tex, posx, posy, posz);
        float4 grad = tex3D(gradientsTex, posx, posy, posz);
        float4 col = tex1D(transferTex, (sample-c_kparams.transferOffset)*c_kparams.transferScale);
        col.w *= c_kparams.density;

        //decrease the gradient magnitude some, to avoid
                //running out of alpha too early.
          //      grad.w *= 0.01;

        /*Boundry Enhancement equation*/
        col.w *= compute_boundary(grad.w);
        
        // pre-multiply alpha
        col.x *= col.w;
        col.y *= col.w;
        col.z *= col.w;

        sum = sum + col*(1.0f - sum.w); // "over" operator for front-to-back blending

        // exit early if opaque
        if (sum.w > opacityThreshold)
            break;

        t += tstep;
        if (t > tfar) break;

        pos += step;
    }
    sum *= c_kparams.brightness;

    // write output color
    d_output[y*c_kparams.imageW + x] = rgbaFloatToInt(sum);
}
/* Distance Color Blending
 * Technique Adpated from: Volume Illustration: Non-Photorealistic Rendering
 *                         of Volume Models by Ebert and Rheingans
 * Makes colors which are farther from the viewer darker, they recede in the distance
 */
__global__ void
d_render_dcb(uint *d_output)
{
    const int maxSteps = 500;
    const float tstep = 0.01f;
    const float opacityThreshold = 0.95f;
    const float3 boxMin = make_float3(-1.0f, -1.0f, -1.0f);
    const float3 boxMax = make_float3(1.0f, 1.0f, 1.0f);

    uint x = blockIdx.x*blockDim.x + threadIdx.x;
    uint y = blockIdx.y*blockDim.y + threadIdx.y;
    if ((x >= c_kparams.imageW) || (y >= c_kparams.imageH)) return;

    float u = (x / (float) c_kparams.imageW)*2.0f-1.0f;
    float v = (y / (float) c_kparams.imageH)*2.0f-1.0f;

    // calculate eye ray in world space
    Ray eyeRay;
    eyeRay.o = make_float3(mul(c_invViewMatrix, make_float4(0.0f, 0.0f, 0.0f, 1.0f)));
    eyeRay.d = normalize(make_float3(u, v, -2.0f));
    eyeRay.d = mul(c_invViewMatrix, eyeRay.d);

    // find intersection with box
    float tnear, tfar;
    int hit = intersectBox(eyeRay, boxMin, boxMax, &tnear, &tfar);
    if (!hit) return;
    if (tnear < 0.0f) tnear = 0.0f;     // clamp to near plane

    // march along ray from front to back, accumulating color
    float4 sum = make_float4(0.0f);
    float t = tnear;
    float3 pos = eyeRay.o + eyeRay.d*tnear;
    float3 step = eyeRay.d*tstep;
   
    float4 cb = make_float4(0.f, 0.f, 0.15f, 1.f);
   

    for(int i=0; i<maxSteps; i++) {
        // read from 3D texture
        // remap position to [0, 1] coordinates
        float posx=  pos.x*0.5f+0.5f;
        float posy=  pos.y*0.5f+0.5f;
        float posz= pos.z*0.5f+0.5f;
        float sample = tex3D(tex, posx, posy, posz);

        // lookup in transfer function texture
        float4 col = tex1D(transferTex, (sample-c_kparams.transferOffset)*c_kparams.transferScale);
        col.w *= c_kparams.density;
        //calculate the fractional distance through the volume
       float dv=(i/(float)maxSteps);
	   //raise the distance by kde
        dv = pow(dv, c_kparams.kde);

	    //modify the sample color to darken with distance
        col.x=(1-c_kparams.kds*dv)*col.x+(c_kparams.kds*dv*cb.x);
        col.y=(1-c_kparams.kds*dv)*col.y+(c_kparams.kds*dv*cb.y);
        col.z=(1-c_kparams.kds*dv)*col.z+(c_kparams.kds*dv*cb.z);

        // pre-multiply alpha
        col.x *= col.w;
        col.y *= col.w;
        col.z *= col.w;

        // "over" operator for front-to-back blending
        sum = sum + col*(1.0f - sum.w);

        // exit early if opaque
        if (sum.w > opacityThreshold)
            break;

        t += tstep;
        if (t > tfar) break;

        pos += step;
    }
    sum *= c_kparams.brightness;

    // write output color
    d_output[y*c_kparams.imageW + x] = rgbaFloatToInt(sum);
}
/* MIDA
 * Technique Adpated from: Instant Volume Visualization
 * Modulates between DVR and MIP. Accumulates in a way to resemble MIP
 */
__global__ void
d_render_mida(uint *d_output)
{
    const int maxSteps = 500;
    const float tstep = 0.01f;
//    const float opacityThreshold = 0.95f;
    const float3 boxMin = make_float3(-1.0f, -1.0f, -1.0f);
    const float3 boxMax = make_float3(1.0f, 1.0f, 1.0f);

    uint x = blockIdx.x*blockDim.x + threadIdx.x;
    uint y = blockIdx.y*blockDim.y + threadIdx.y;
    if ((x >= c_kparams.imageW) || (y >= c_kparams.imageH)) return;

    float u = (x / (float) c_kparams.imageW)*2.0f-1.0f;
    float v = (y / (float) c_kparams.imageH)*2.0f-1.0f;

    // calculate eye ray in world space
    Ray eyeRay;
    eyeRay.o = make_float3(mul(c_invViewMatrix, make_float4(0.0f, 0.0f, 0.0f, 1.0f)));
    eyeRay.d = normalize(make_float3(u, v, -2.0f));
    eyeRay.d = mul(c_invViewMatrix, eyeRay.d);

    // find intersection with box
    float tnear, tfar;
    int hit = intersectBox(eyeRay, boxMin, boxMax, &tnear, &tfar);
    if (!hit) return;
    if (tnear < 0.0f) tnear = 0.0f;     // clamp to near plane

    float4 sum = make_float4(0.0f);
    float t = tnear;
    float3 pos = eyeRay.o + eyeRay.d*tnear;
    float3 step = eyeRay.d*tstep;

    float4 max_rgb = make_float4(0.f);
    float delta=0;
    float mida_max=0;
    float beta=0;
    
    for(int i=0; i<maxSteps; i++) {

        float posx = pos.x*0.5f+0.5f;
        float posy = pos.y*0.5f+0.5f;
        float posz = pos.z*0.5f+0.5f;

        float sample = tex3D(tex, posx, posy, posz);
        float4 col   = tex1D(transferTex, (sample-c_kparams.transferOffset)*c_kparams.transferScale);

        /*Calculate delta, if a new maximum is encountered, save it
	 *as the new maximum (in mida_max)
	 * and save the color of the maximum sample. If not a max, delta=0
	 */
        if(sample>mida_max)
        {
            delta=sample-mida_max;
            mida_max=sample;
	    max_rgb=col;
        }
        else
        {
            delta=0;
        }
    	col.w *= c_kparams.density;
     	/*If the gamma is less than 0 will have a combo of  MIDA and
	 *DVR forms. When gamma is 0 it will be fully DVR
	 */
        if(c_kparams.mida_gamma<0)
        {	
            beta=1-(delta*(1+c_kparams.mida_gamma));
        }

        else
        {
            beta=1-delta;
        }
	//accumulate the color and the alpha value by interpolating with beta as a weight
        sum.x=beta*(sum.x)+(1.0f-beta*sum.w)*col.w*col.x;
        sum.y=beta*(sum.y)+(1.0f-beta*sum.w)*col.w*col.y;
        sum.z=beta*(sum.z)+(1.0f-beta*sum.w)*col.w*col.z;
        sum.w=beta*sum.w+(1-beta*sum.w)*col.w;

        t += tstep;
        if (t > tfar) break;

        pos += step;
    }
     /*If the gamma is greater than 0 we must interpolate between the MIDA and
	*MIP forms. When gamma is 1 it will be fully MIP*/
     if(c_kparams.mida_gamma >0)
     {
	    sum.x=(1.0-c_kparams.mida_gamma)*(sum.x*sum.w)+(c_kparams.mida_gamma)*(max_rgb.x*max_rgb.w);
	    sum.y=(1.0-c_kparams.mida_gamma)*(sum.y*sum.w)+(c_kparams.mida_gamma)*(max_rgb.y*max_rgb.w);
	    sum.z=(1.0-c_kparams.mida_gamma)*(sum.z*sum.w)+(c_kparams.mida_gamma)*(max_rgb.z*max_rgb.w);
	    sum.w=(1.0-c_kparams.mida_gamma)*(sum.w*sum.w)+(c_kparams.mida_gamma)*(max_rgb.w);
     }

    // write output color
    d_output[y*c_kparams.imageW + x] = rgbaFloatToInt(sum);
}

/* Tone Shading
 * Technique Adpated from: Volume Illustration: Non-Photorealistic Rendering
 *                         of Volume Models by Ebert and Rheingans
 * Volume appears to be lit from a warm light, objects facing the light are yellow
 * 	  objects away from the light are blue
 */
__global__ void
d_render_tone(uint *d_output)
{
    const int maxSteps = 500;
    const float tstep = 0.01f;
    const float opacityThreshold = 0.95f;
    const float3 boxMin = make_float3(-1.0f, -1.0f, -1.0f);
    const float3 boxMax = make_float3(1.0f, 1.0f, 1.0f);

    uint x = blockIdx.x*blockDim.x + threadIdx.x;
    uint y = blockIdx.y*blockDim.y + threadIdx.y;
    if ((x >= c_kparams.imageW) || (y >= c_kparams.imageH)) return;

    float u = (x / (float) c_kparams.imageW)*2.0f-1.0f;
    float v = (y / (float) c_kparams.imageH)*2.0f-1.0f;

    // calculate eye ray in world space
    Ray eyeRay;
    eyeRay.o = make_float3(mul(c_invViewMatrix, make_float4(0.0f, 0.0f, 0.0f, 1.0f)));
    eyeRay.d = normalize(make_float3(u, v, -2.0f));
    eyeRay.d = mul(c_invViewMatrix, eyeRay.d);

    // find intersection with box
    float tnear, tfar;
    int hit = intersectBox(eyeRay, boxMin, boxMax, &tnear, &tfar);
    if (!hit) return;
    if (tnear < 0.0f) tnear = 0.0f;     // clamp to near plane

    float4 sum = make_float4(0.0f);
    float t = tnear;
    float3 pos = eyeRay.o + eyeRay.d*tnear;
    float3 step = eyeRay.d*tstep;

    for(int i=0; i<maxSteps; i++) {

        float posx=  pos.x*0.5f+0.5f;
        float posy=  pos.y*0.5f+0.5f;
        float posz=  pos.z*0.5f+0.5f;

        float sample = tex3D(tex, posx, posy, posz);
        float4 col   = tex1D(transferTex, (sample-c_kparams.transferOffset)*c_kparams.transferScale);
        float4 grad =  tex3D(gradientsTex, posx, posy, posz);
        col.w *= c_kparams.density;
        float4 col_tone = compute_tone(col,grad,pos);
        col.x = col_tone.x;
        col.y = col_tone.y;
        col.z = col_tone.z;

        //pre-multiply alpha
        col.x *= col.w;
        col.y *= col.w;
        col.z *= col.w;

        sum = sum + col*(1.0f - sum.w);  // "over" operator for front-to-back blending

        if (sum.w > opacityThreshold)
        	break;

        t += tstep;
        if (t > tfar) break;

        pos += step;
    }
    sum *= c_kparams.brightness;

    // write output color
    d_output[y*c_kparams.imageW + x] = rgbaFloatToInt(sum);
}

__global__ void
d_render_tone_class(uint *d_output)
{
    const int maxSteps = 500;
    const float tstep = 0.01f;
    const float opacityThreshold = 0.95f;
    const float3 boxMin = make_float3(-1.0f, -1.0f, -1.0f);
    const float3 boxMax = make_float3(1.0f, 1.0f, 1.0f);

    uint x = blockIdx.x*blockDim.x + threadIdx.x;
    uint y = blockIdx.y*blockDim.y + threadIdx.y;
    if ((x >= c_kparams.imageW) || (y >= c_kparams.imageH)) return;

    float u = (x / (float) c_kparams.imageW)*2.0f-1.0f;
    float v = (y / (float) c_kparams.imageH)*2.0f-1.0f;

    // calculate eye ray in world space
    Ray eyeRay;
    eyeRay.o = make_float3(mul(c_invViewMatrix, make_float4(0.0f, 0.0f, 0.0f, 1.0f)));
    eyeRay.d = normalize(make_float3(u, v, -2.0f));
    eyeRay.d = mul(c_invViewMatrix, eyeRay.d);

    // find intersection with box
    float tnear, tfar;
    int hit = intersectBox(eyeRay, boxMin, boxMax, &tnear, &tfar);
    if (!hit) return;
    if (tnear < 0.0f) tnear = 0.0f;     // clamp to near plane

    float4 sum = make_float4(0.0f);
    float t = tnear;
    float3 pos = eyeRay.o + eyeRay.d*tnear;
    float3 step = eyeRay.d*tstep;

    for(int i=0; i<maxSteps; i++) {

        float posx=  pos.x*0.5f+0.5f;
        float posy=  pos.y*0.5f+0.5f;
        float posz=  pos.z*0.5f+0.5f;

        float sample = tex3D(tex, posx, posy, posz);
        float4 col   = tex1D(transferTex, (sample-c_kparams.transferOffset)*c_kparams.transferScale);
        float4 gradvec =  tex3D(gradientsTex, posx, posy, posz);

        col.w *= c_kparams.density;
        float alpha = classify(sample, gradvec.w);
        if (alpha==0){
            t += tstep;
            if (t > tfar) break;

            pos += step;
            continue;
        }

        float4 col_tone = compute_tone(col,gradvec,pos);
        
        col.x = col_tone.x;
        col.y = col_tone.y;
        col.z = col_tone.z;
        col.w *= alpha * compute_boundary(gradvec.w);


        //pre-multiply alpha
        col.x *= col.w;
        col.y *= col.w;
        col.z *= col.w;

        sum = sum + col*(1.0f - sum.w);  // "over" operator for front-to-back blending

        if (sum.w > opacityThreshold)
            break;

        t += tstep;
        if (t > tfar) break;

        pos += step;
    }
    sum *= c_kparams.brightness;

    // write output color
    d_output[y*c_kparams.imageW + x] = rgbaFloatToInt(sum);
}

__global__ void
d_render_dcb_lighting_class(uint *d_output)
{
    const int maxSteps = 500;
    const float tstep = 0.01f;
    const float opacityThreshold = 0.95f;
    const float3 boxMin = make_float3(-1.0f, -1.0f, -1.0f);
    const float3 boxMax = make_float3(1.0f, 1.0f, 1.0f);

    uint x = blockIdx.x*blockDim.x + threadIdx.x;
    uint y = blockIdx.y*blockDim.y + threadIdx.y;
    if ((x >= c_kparams.imageW) || (y >= c_kparams.imageH)) return;

    float u = (x / (float) c_kparams.imageW)*2.0f-1.0f;
    float v = (y / (float) c_kparams.imageH)*2.0f-1.0f;

    // calculate eye ray in world space
    Ray eyeRay;
    eyeRay.o = make_float3(mul(c_invViewMatrix, make_float4(0.0f, 0.0f, 0.0f, 1.0f)));
    eyeRay.d = normalize(make_float3(u, v, -2.0f));
    eyeRay.d = mul(c_invViewMatrix, eyeRay.d);

    // find intersection with box
    float tnear, tfar;
    int hit = intersectBox(eyeRay, boxMin, boxMax, &tnear, &tfar);
    if (!hit) return;
    if (tnear < 0.0f) tnear = 0.0f;     // clamp to near plane

    // march along ray from front to back, accumulating color
    float4 sum = make_float4(0.0f);
    float t = tnear;
    float3 pos = eyeRay.o + eyeRay.d*tnear;
    float3 step = eyeRay.d*tstep;

    float4 cb = make_float4(0.f, 0.f, 0.15f, 1.f);
    float3 lightPos = make_float3(c_kparams.lightx, c_kparams.lighty, c_kparams.lightz);

    for(int i=0; i<maxSteps; i++) {
        // read from 3D texture
        // remap position to [0, 1] coordinates
        float posx=  pos.x*0.5f+0.5f;
        float posy=  pos.y*0.5f+0.5f;
        float posz= pos.z*0.5f+0.5f;
        float sample = tex3D(tex, posx, posy, posz);

        // lookup in transfer function texture
        float4 col = tex1D(transferTex, (sample-c_kparams.transferOffset)*c_kparams.transferScale);
        float4 gradvec = tex3D(gradientsTex, posx, posy, posz);
        float3 grad3 = make_float3(gradvec.x, gradvec.y, gradvec.z);
        col.w *= c_kparams.density;
        
        //calculate the fractional distance through the volume
        float dv=(i/(float)maxSteps);
        dv = pow(dv, c_kparams.kde);

        //modify the sample color to darken with distance
        col.x=(1-c_kparams.kds*dv)*col.x+(c_kparams.kds*dv*cb.x);
        col.y=(1-c_kparams.kds*dv)*col.y+(c_kparams.kds*dv*cb.y);
        col.z=(1-c_kparams.kds*dv)*col.z+(c_kparams.kds*dv*cb.z);

        float3 light = compute_light(grad3, normalize(lightPos-pos), eyeRay.d);        
        float alpha = classify(sample, gradvec.w);
        if (alpha==0){
            t += tstep;
            if (t > tfar) break;

            pos += step;
            continue;
        }
        col.w *= alpha;
        col.x *= light.x;
        col.y *= light.y;
        col.z *= light.z;


        // pre-multiply alpha
        col.x *= col.w;
        col.y *= col.w;
        col.z *= col.w;

        // "over" operator for front-to-back blending
        sum = sum + col*(1.0f - sum.w);

        // exit early if opaque
        if (sum.w > opacityThreshold)
            break;

        t += tstep;
        if (t > tfar) break;

        pos += step;
    }
    sum *= c_kparams.brightness;

    // write output color
    d_output[y*c_kparams.imageW + x] = rgbaFloatToInt(sum);
}

__global__ void
d_render_depth_buffer(float *d_output, float *d_dbuf)
{
    const int maxSteps = 500;
    const float tstep = 0.01f;
    const float opacityThreshold = 0.95f;
    const float3 boxMin = make_float3(-1.0f, -1.0f, -1.0f);
    const float3 boxMax = make_float3(1.0f, 1.0f, 1.0f);

    uint x = blockIdx.x*blockDim.x + threadIdx.x;
    uint y = blockIdx.y*blockDim.y + threadIdx.y;
    if ((x >= c_kparams.imageW) || (y >= c_kparams.imageH)) return;

    float u = (x / (float) c_kparams.imageW)*2.0f-1.0f;
    float v = (y / (float) c_kparams.imageH)*2.0f-1.0f;

    // calculate eye ray in world space
    Ray eyeRay;
    eyeRay.o = make_float3(mul(c_invViewMatrix, make_float4(0.0f, 0.0f, 0.0f, 1.0f)));
    eyeRay.d = normalize(make_float3(u, v, -2.0f));
    eyeRay.d = mul(c_invViewMatrix, eyeRay.d);

    // find intersection with box
    float tnear, tfar;
    int hit = intersectBox(eyeRay, boxMin, boxMax, &tnear, &tfar);
    if (!hit) return;
    if (tnear < 0.0f) tnear = 0.0f;     // clamp to near plane

    // march along ray from front to back, accumulating color
    float4 sum = make_float4(0.0f);
    float t = tnear;
    float3 pos = eyeRay.o + eyeRay.d*tnear;
    float3 step = eyeRay.d*tstep;

    for(int i=0; i<maxSteps; i++) {
        float posx=  pos.x*0.5f+0.5f;
        float posy=  pos.y*0.5f+0.5f;
        float posz=  pos.z*0.5f+0.5f;

        float sample = tex3D(tex, posx, posy, posz);
        float4 col = tex1D(transferTex, (sample-c_kparams.transferOffset)*c_kparams.transferScale);
        col.w *= c_kparams.density;

        // pre-multiply alpha
        col.x *= col.w;
        col.y *= col.w;
        col.z *= col.w;

        sum = sum + col*(1.0f - sum.w); // "over" operator for front-to-back blending

        // exit early if opaque
        if (sum.w > opacityThreshold)
            break;

        t += tstep;
        if (t > tfar) break;

        pos += step;
    }
    sum *= c_kparams.brightness;

    // write output color
    d_output[y*c_kparams.imageW + x] = rgbaFloatToInt(sum);
}


extern "C"
void render_kernel(dim3 gridSize, dim3 blockSize, 
                   uint *d_output, float *d_dbuf, 
                   int method)
{
    switch (method){
    case 1:
        d_render_orig<<<gridSize, blockSize>>>(d_output);
        break;
    case 2:
    	d_render_classification<<<gridSize, blockSize>>>(d_output);
    	break;
    case 3:
        d_render_bound<<<gridSize, blockSize>>>(d_output);
        break;
    case 4:
        d_render_dcb<<<gridSize, blockSize>>>(d_output);
        break;
    case 5:
        d_render_tone<<<gridSize, blockSize>>>(d_output);
        break;
    case 6:
        d_render_mida<<<gridSize, blockSize>>>(d_output);
        break;
    case 7:
        d_render_sil<<<gridSize,blockSize>>>(d_output);
        break;
    case 8:
        d_render_tone_class<<<gridSize,blockSize>>>(d_output);
        break;
    case 9:
        d_render_dcb_lighting_class<<<gridSize,blockSize>>>(d_output);
        break;
//     case 7:
//         d_render_depth_buffer<<<gridSize,blockSize>>>(d_depthBuffer);
//         break;
//     case 100:
//         d_compute_depth_buffer<<<gridSize, blockSize>>>(d_depthBuffer);
//         break;
    default:
       	d_render_orig<<<gridSize, blockSize>>>(d_output);
       	break;
    }
}
#endif // #ifndef _VOLUMERENDER_KERNEL_CU_



